SummaryPlot

This R script is used for generating summary plot for the Geneconv project

1, Read in tables

rm(list=ls())  # clean up workspace
path.HKY <- "/Users/xji3/Genconv/NewClusterPackRun/NewPackageNewRun/"
path.MG94 <- "/Users/xji3/FromCluster03192015/NewPackageNewRun/"
path.HKY <- path.MG94
#path.MG94 = path.HKY
#path <- "G:/Geneconv/NewClusterPackRun/NewPackageNewRun/"
#HKY_clock_summary <- "HKY_clock_summary"
summary.list.MG94 <- c( 
                  "MG94_clock_summary",
                  "MG94_nonclock_summary", 
                  "Force_MG94_clock_summary", 
                  "Force_MG94_nonclock_summary"
                  )
summary.list.HKY <- c( "HKY_nonclock_summary",
                  "HKY_clock_summary",
                  "Force_HKY_clock_summary", 
                 "Force_HKY_nonclock_summary"
                  )
#summary.list <- c("Force_MG94_clock_summary")
for (target.summary in summary.list.MG94){
  summary_file <- paste(path.MG94, target.summary, '.txt', sep = '')
  all <- readLines(summary_file, n = -1)
  col.names <- strsplit(all[1], ' ')[[1]][-1]
  row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
  summary_mat <- as.matrix(read.table(summary_file, 
                                        row.names = row.names, 
                                        col.names = col.names))
  assign(target.summary, summary_mat)
}

for (target.summary in summary.list.HKY){
  summary_file <- paste(path.HKY, target.summary, '.txt', sep = '')
  all <- readLines(summary_file, n = -1)
  col.names <- strsplit(all[1], ' ')[[1]][-1]
  row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
  summary_mat <- as.matrix(read.table(summary_file, 
                                        row.names = row.names, 
                                        col.names = col.names))
  assign(target.summary, summary_mat)
}

# read in directional model summary
summary.list.Dir <- c( 
                  "MG94_clock_summary",
                  "MG94_nonclock_summary",
                  "HKY_nonclock_summary",
                  "HKY_clock_summary"
                  )
for (target.summary in summary.list.Dir){
  summary_file <- paste(path.MG94, "Dir_" , target.summary, '.txt', sep = '')
  all <- readLines(summary_file, n = -1)
  col.names <- strsplit(all[1], ' ')[[1]][-1]
  row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
  summary_mat <- as.matrix(read.table(summary_file, 
                                        row.names = row.names, 
                                        col.names = col.names))
  assign(paste("Dir_", target.summary, sep = ''), summary_mat)
  write.table(gsub("_", ".", paste("Dir_", target.summary, sep = '')), 
              paste(path.MG94, "Dir_", target.summary, sep = ''))
}

# read in gBGC model summary
summary.list.gBGC <- c( 
                  "MG94_clock_summary",
                  "MG94_nonclock_summary",
                  "HKY_nonclock_summary",
                  "HKY_clock_summary"
                  )
for (target.summary in summary.list.gBGC){
  summary_file <- paste(path.MG94, "gBGC_" , target.summary, '.txt', sep = '')
  all <- readLines(summary_file, n = -1)
  col.names <- strsplit(all[1], ' ')[[1]][-1]
  row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
  summary_mat <- as.matrix(read.table(summary_file, 
                                        row.names = row.names, 
                                        col.names = col.names))
  assign(paste("gBGC_", target.summary, sep = ''), summary_mat)
  write.table(gsub("_", ".", paste("gBGC_", target.summary, sep = '')), 
              paste(path.MG94, "gBGC_", target.summary, sep = ''))
}

# read in Dir gBGC model summary
summary.list.gBGC <- c( 
                  "MG94_clock_summary",
                  "MG94_nonclock_summary",
                  "HKY_nonclock_summary",
                  "HKY_clock_summary"
                  )
for (target.summary in summary.list.gBGC){
  summary_file <- paste(path.MG94, "Dir_gBGC_" , target.summary, '.txt', sep = '')
  all <- readLines(summary_file, n = -1)
  col.names <- strsplit(all[1], ' ')[[1]][-1]
  row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
  summary_mat <- as.matrix(read.table(summary_file, 
                                        row.names = row.names, 
                                        col.names = col.names))
  assign(paste("Dir_gBGC_", target.summary, sep = ''), summary_mat)
  write.table(gsub("_", ".", paste("Dir_gBGC_", target.summary, sep = '')), 
              paste(path.MG94, "Dir_gBGC_", target.summary, sep = ''))
}

Now generate summary file of only pairs that have all cases finished in HKY or MG94 models.

# HKY
HKY.pair.names <- intersect(intersect(intersect(colnames(HKY_clock_summary), 
                            colnames(HKY_nonclock_summary)),
                            colnames(Force_HKY_clock_summary)),
                            colnames(Force_HKY_nonclock_summary))
HKY.clock.filtered <- HKY_clock_summary[, HKY.pair.names]
HKY.nonclock.filtered <- HKY_nonclock_summary[, HKY.pair.names]
HKY.Force.clock.filtered <- Force_HKY_clock_summary[, HKY.pair.names]
HKY.Force.nonclock.filtered <- Force_HKY_nonclock_summary[, HKY.pair.names]

write.table(HKY.clock.filtered, paste( path.HKY, "HKY_clock_filtered", sep = ""))
write.table(HKY.nonclock.filtered, paste( path.HKY, "HKY_nonclock_filtered", sep = ""))
write.table(HKY.Force.clock.filtered, paste( path.HKY, "HKY_Force_clock_filtered", sep = ""))
write.table(HKY.Force.nonclock.filtered, paste( path.HKY, "HKY_Force_nonclock_filtered", sep = ""))



# MG94
MG94.pair.names <- intersect(intersect(intersect(colnames(MG94_clock_summary), 
                            colnames(MG94_nonclock_summary)),
                            colnames(Force_MG94_clock_summary)),
                            colnames(Force_MG94_nonclock_summary))
MG94.clock.filtered <- MG94_clock_summary[, MG94.pair.names]
MG94.nonclock.filtered <- MG94_nonclock_summary[, MG94.pair.names]
MG94.Force.clock.filtered <- Force_MG94_clock_summary[, MG94.pair.names]
MG94.Force.nonclock.filtered <- Force_MG94_nonclock_summary[, MG94.pair.names]

write.table(MG94.clock.filtered, paste( path.MG94, "MG94_clock_filtered", sep = ""))
write.table(MG94.nonclock.filtered, paste( path.MG94, "MG94_nonclock_filtered", sep = ""))
write.table(MG94.Force.clock.filtered, paste( path.MG94, "MG94_Force_clock_filtered", sep = ""))
write.table(MG94.Force.nonclock.filtered, paste( path.MG94, "MG94_Force_nonclock_filtered", sep = ""))

Now analyze the results

First, show the loglikelihood improvement for each model with/without tau

# HKY nonclock
(HKY.nonclock.filtered - HKY.Force.nonclock.filtered)[2,]

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##           45.32           54.37          211.72          135.89 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##          152.26          165.28           42.03           46.49 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##           90.04           42.93           60.11           73.59 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##          156.46           44.68           13.36           94.14 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##          230.93           89.89           27.77           74.34 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##          120.41          155.39           96.07          158.20 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##           96.69          109.92           16.80          114.89 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##         1336.94           94.99           85.48          108.80 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##          203.71          232.13           22.94           17.39

# HKY clock
(HKY.clock.filtered - HKY.Force.clock.filtered)[2,]

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##          45.671          47.253         169.574         135.283 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##         149.475         145.207          41.107          39.588 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##          77.690          34.149          40.567          63.727 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##         152.259          32.103          15.203          94.121 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##         231.875          92.031          11.503          74.361 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##         119.181         163.761          93.715         121.540 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##          75.862          92.364          17.778         113.481 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##        1262.228          76.723          66.710         118.843 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##         201.896         229.128          26.628           7.088

# # MG94 nonclock
(MG94.nonclock.filtered - MG94.Force.nonclock.filtered)[2,]

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##          16.983           6.754           5.705          56.476 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##          75.763          22.251          16.404          11.801 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##          58.607          33.775           5.514           3.031 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##          14.255           8.070           3.765          37.228 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##          54.544          36.620           5.434          21.596 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##          56.946         109.590          32.040          52.565 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##          18.334          30.024           0.892          33.229 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##         282.233          30.934           7.521          83.647 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##         139.211         130.699          14.646           3.204

# # MG94 clock
(MG94.clock.filtered - MG94.Force.clock.filtered)[2,]

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##         17.4243         15.5316         39.4858         68.2596 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##         79.4879         20.3225         16.7182         11.3831 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##         53.3016         30.6935          3.1882          3.4542 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##         15.7566          6.6699          6.0741         39.0434 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##         57.5931         38.7124          1.8099         23.4107 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##         66.8919        118.0286         30.5241         48.0554 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##         13.5553         28.8624          0.8200         49.5502 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##        480.0589         25.3489          5.3776         92.3178 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##        146.4337        140.9595         18.1782          0.6007

Then, show loglikelihood improvement for other nested models

# HKY nonclock & clock
(HKY.nonclock.filtered[2, ] - HKY.clock.filtered[2,])

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##           5.274          14.824          66.101           7.517 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##           8.969          53.416          10.808          36.218 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##          18.948          56.902          61.573          30.532 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##          31.128          38.191          44.197          12.081 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##          30.547           7.267         135.687          12.446 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##          15.795           5.421          14.239         115.644 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##          86.961          33.687          29.250           8.811 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##         119.681          65.794          54.428           4.112 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##          10.140          10.057          52.538          64.495

# MG94 nonclock & clock
(MG94.nonclock.filtered[2, ] - MG94.clock.filtered[2,])

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##           4.037          -3.046         -11.566          -5.696 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##           3.010          23.198           6.948          17.283 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##          11.254          33.909          31.542          13.044 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##          13.385          14.248          43.391           8.979 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##          13.242           5.237          97.800           6.140 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##           2.383           5.864          12.064          74.575 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##          50.585          21.955          18.096          -9.918 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##        -138.441          34.134          23.873           3.269 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##           1.445          -2.976          32.539          49.944

Then, show the edge specific tau estimates (posterior expected number of geneconv events / posterior expected time in heterogeneous states of each branch)

HKY.nonclock.filtered[21:32, ]

##                       YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W
## (N0,N1,tau)                     8.017           7.160           2.942
## (N0,kluyveri,tau)               0.000           0.000           0.000
## (N1,N2,tau)                    10.018           8.120           3.116
## (N1,castellii,tau)              9.960           7.194           3.169
## (N2,N3,tau)                     5.498           8.109           2.829
## (N2,bayanus,tau)                8.304           8.054           3.168
## (N3,N4,tau)                     5.516           6.566           3.245
## (N3,kudriavzevii,tau)           4.163           7.662           3.551
## (N4,N5,tau)                    12.566           6.636           3.243
## (N4,mikatae,tau)                7.693           8.279           3.045
## (N5,cerevisiae,tau)             6.501          11.617           3.404
## (N5,paradoxus,tau)              8.500           4.073           2.922
##                       YML026C_YDR450W YNL301C_YOL120C YHR106W_YDR353W
## (N0,N1,tau)                     14.08          10.597           1.801
## (N0,kluyveri,tau)                0.00           0.000           0.000
## (N1,N2,tau)                     13.74          10.203           1.573
## (N1,castellii,tau)              18.48          11.180           2.227
## (N2,N3,tau)                     15.18          11.790           1.652
## (N2,bayanus,tau)                17.00          14.142           1.820
## (N3,N4,tau)                      6.15           8.830           1.743
## (N3,kudriavzevii,tau)           15.50          10.152           1.682
## (N4,N5,tau)                     17.23          11.714           1.677
## (N4,mikatae,tau)                13.87           7.845           1.791
## (N5,cerevisiae,tau)             19.41          17.769           2.122
## (N5,paradoxus,tau)              16.32           6.928           1.706
##                       YIL057C_YER067W YPL087W_YBR183W YNL069C_YIL133C
## (N0,N1,tau)                    0.9864          0.4711           5.100
## (N0,kluyveri,tau)              0.0000          0.0000           0.000
## (N1,N2,tau)                    0.9380          0.4248           5.061
## (N1,castellii,tau)             1.1657          0.5068           6.308
## (N2,N3,tau)                    0.9199          0.3093           6.175
## (N2,bayanus,tau)               0.7904          0.4962           5.878
## (N3,N4,tau)                    0.7103          0.4493           6.438
## (N3,kudriavzevii,tau)          0.9511          0.4387           3.809
## (N4,N5,tau)                    0.6357          0.4751           4.839
## (N4,mikatae,tau)               1.1072          0.4835           3.038
## (N5,cerevisiae,tau)            0.7711          0.5487           4.259
## (N5,paradoxus,tau)             1.1634          0.4167           4.386
##                       YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W
## (N0,N1,tau)                    0.2834          0.7241          0.5580
## (N0,kluyveri,tau)              0.0000          0.0000          0.0000
## (N1,N2,tau)                    0.2611          0.6626          0.5162
## (N1,castellii,tau)             0.3080          0.7194          0.5426
## (N2,N3,tau)                    0.2119          0.6608          0.4071
## (N2,bayanus,tau)               0.2641          0.7238          0.5770
## (N3,N4,tau)                    0.3089          0.7189          0.6951
## (N3,kudriavzevii,tau)          0.3054          0.6704          0.6129
## (N4,N5,tau)                    0.2915          0.6790          0.6485
## (N4,mikatae,tau)               0.2731          0.8823          0.5481
## (N5,cerevisiae,tau)            0.3072          0.8881          0.6488
## (N5,paradoxus,tau)             0.2865          0.7330          0.5631
##                       YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W
## (N0,N1,tau)                     3.854          0.4731          0.3706
## (N0,kluyveri,tau)               0.000          0.0000          0.0000
## (N1,N2,tau)                     3.265          0.4647          0.3192
## (N1,castellii,tau)              4.823          0.4979          0.5583
## (N2,N3,tau)                     3.452          0.4399          0.3054
## (N2,bayanus,tau)                3.335          0.3921          0.3522
## (N3,N4,tau)                     3.519          0.6365          0.4992
## (N3,kudriavzevii,tau)           4.400          0.4855          0.3867
## (N4,N5,tau)                     3.196          0.3820          0.4002
## (N4,mikatae,tau)                4.733          0.4949          0.4117
## (N5,cerevisiae,tau)             3.504          0.4748          0.3476
## (N5,paradoxus,tau)              3.695          0.4279          0.3069
##                       YMR143W_YDL083C YPR159W_YGR143W YJL177W_YKL180W
## (N0,N1,tau)                    12.341          1.1321          10.082
## (N0,kluyveri,tau)               0.000          0.0000           0.000
## (N1,N2,tau)                    12.444          1.1459           9.838
## (N1,castellii,tau)             14.628          1.1930          10.315
## (N2,N3,tau)                     5.070          0.9358           7.493
## (N2,bayanus,tau)               20.536          1.1858          14.605
## (N3,N4,tau)                    15.679          0.9418          12.588
## (N3,kudriavzevii,tau)           9.011          1.0776           6.409
## (N4,N5,tau)                    13.233          1.1149          15.819
## (N4,mikatae,tau)               11.703          1.1808           9.909
## (N5,cerevisiae,tau)            16.503          0.9481           9.672
## (N5,paradoxus,tau)             17.976          1.2592          10.649
##                       YGL133W_YPL216W YJR048W_YEL039C YBR191W_YPL079W
## (N0,N1,tau)                    0.1251           2.315           14.46
## (N0,kluyveri,tau)              0.0000           0.000            0.00
## (N1,N2,tau)                    0.1279           2.084           14.12
## (N1,castellii,tau)             0.1320           2.817           15.50
## (N2,N3,tau)                    0.1205           1.868           16.26
## (N2,bayanus,tau)               0.1054           2.283           20.45
## (N3,N4,tau)                    0.1243           2.401           17.97
## (N3,kudriavzevii,tau)          0.1233           2.067           12.28
## (N4,N5,tau)                    0.1327           2.101           15.60
## (N4,mikatae,tau)               0.1223           2.059           20.27
## (N5,cerevisiae,tau)            0.1247           2.426           11.86
## (N5,paradoxus,tau)             0.1079           2.430           20.01
##                       YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C
## (N0,N1,tau)                     23.52           7.993          0.6125
## (N0,kluyveri,tau)                0.00           0.000          0.0000
## (N1,N2,tau)                     26.13           7.428          0.6319
## (N1,castellii,tau)              20.18           7.910          0.6122
## (N2,N3,tau)                     27.82           7.715          0.5484
## (N2,bayanus,tau)                37.03          10.498          0.6391
## (N3,N4,tau)                     36.99           6.876          0.5646
## (N3,kudriavzevii,tau)           15.12           7.088          0.5623
## (N4,N5,tau)                     39.47           9.724          0.5540
## (N4,mikatae,tau)                32.30           8.065          0.6168
## (N5,cerevisiae,tau)             27.16          10.640          0.6068
## (N5,paradoxus,tau)              51.57           8.115          0.6763
##                       YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C
## (N0,N1,tau)                    0.6618           1.190          0.2080
## (N0,kluyveri,tau)              0.0000           0.000          0.0000
## (N1,N2,tau)                    0.5784           1.230          0.2005
## (N1,castellii,tau)             0.6638           1.207          0.1954
## (N2,N3,tau)                    0.6002           1.040          0.2099
## (N2,bayanus,tau)               0.6426           1.259          0.1939
## (N3,N4,tau)                    0.6717           1.224          0.1926
## (N3,kudriavzevii,tau)          0.6503           1.107          0.2324
## (N4,N5,tau)                    0.6475           1.268          0.2713
## (N4,mikatae,tau)               0.7376           1.021          0.2284
## (N5,cerevisiae,tau)            0.8215           1.231          0.2046
## (N5,paradoxus,tau)             0.6711           1.355          0.2050
##                       YBL087C_YER117W YGL062W_YBR218C YIR033W_YKL020C
## (N0,N1,tau)                    12.021           4.038          0.3405
## (N0,kluyveri,tau)               0.000           0.000          0.0000
## (N1,N2,tau)                    11.927           4.389          0.3492
## (N1,castellii,tau)             15.646           4.039          0.3532
## (N2,N3,tau)                    10.741           4.216          0.3380
## (N2,bayanus,tau)               21.791           5.272          0.3501
## (N3,N4,tau)                    11.055           4.136          0.2873
## (N3,kudriavzevii,tau)          13.864           4.371          0.3302
## (N4,N5,tau)                    18.699           4.520          0.3051
## (N4,mikatae,tau)                8.222           4.272          0.2908
## (N5,cerevisiae,tau)            17.032           5.118          0.3201
## (N5,paradoxus,tau)             12.448           3.986          0.3397
##                       YMR243C_YOR316C YLR333C_YGR027C YMR142C_YDL082W
## (N0,N1,tau)                    0.6858          11.015          15.012
## (N0,kluyveri,tau)              0.0000           0.000           0.000
## (N1,N2,tau)                    0.7021          14.043          14.898
## (N1,castellii,tau)             0.8128          10.546          17.123
## (N2,N3,tau)                    0.5933           8.322          11.765
## (N2,bayanus,tau)               0.6732          20.725          23.207
## (N3,N4,tau)                    0.6194          11.616          16.446
## (N3,kudriavzevii,tau)          0.6137           4.684          18.693
## (N4,N5,tau)                    0.5944           6.069          18.335
## (N4,mikatae,tau)               0.6332           8.513          12.273
## (N5,cerevisiae,tau)            0.6755          10.147           9.302
## (N5,paradoxus,tau)             0.5864           5.497          17.758
##                       YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C
## (N0,N1,tau)                    13.839          0.1455          0.2115
## (N0,kluyveri,tau)               0.000          0.0000          0.0000
## (N1,N2,tau)                    15.687          0.1428          0.2219
## (N1,castellii,tau)             15.345          0.1408          0.2394
## (N2,N3,tau)                    16.285          0.1330          0.1774
## (N2,bayanus,tau)               19.523          0.1336          0.2194
## (N3,N4,tau)                    23.204          0.1435          0.1696
## (N3,kudriavzevii,tau)           9.836          0.1363          0.1661
## (N4,N5,tau)                    32.491          0.1557          0.1883
## (N4,mikatae,tau)               15.836          0.1470          0.1990
## (N5,cerevisiae,tau)            48.908          0.2033          0.1564
## (N5,paradoxus,tau)             23.208          0.1530          0.1919

Now plot Total blen v.s. Tau into different groups (differ by color)

HKY nonclock case

plot(colSums(HKY.Force.nonclock.filtered[9:20, ]), HKY.nonclock.filtered[8, ],
     type = "n", xlab = "Total branch length", ylab = "Tau" )
col.color <- rep("black", dim(HKY.nonclock.filtered)[2])
col.color[(HKY.nonclock.filtered - HKY.Force.nonclock.filtered)[2,] < 50] <- "red"
points(x = colSums(HKY.Force.nonclock.filtered[9:20, ]), y = HKY.nonclock.filtered[8, ],
       type = "p", col = col.color, bg = col.color)
legend("topright", c("lnL improvement > 50", "lnL improvement < 50"),
       lty = c(1, 1), 
       lwd = c(2.5, 2.5),
       col = c("black", "red"))
title("HKY nonclock")

HKY clock case

plot(colSums(HKY.Force.clock.filtered[9:20, ]), HKY.clock.filtered[8, ],
     type = "n", xlab = "Total branch length", ylab = "Tau" )
col.color <- rep("black", dim(HKY.clock.filtered)[2])
col.color[(HKY.clock.filtered - HKY.Force.clock.filtered)[2,] < 50] <- "red"
points(x = colSums(HKY.Force.clock.filtered[9:20, ]), y = HKY.clock.filtered[8, ],
       type = "p", col = col.color, bg = col.color)
legend("topright", c("lnL improvement > 50", "lnL improvement < 50"),
       lty = c(1, 1), 
       lwd = c(2.5, 2.5),
       col = c("black", "red"))
title("HKY clock")

MG94 nonclock case

improvement.lmt <- 20
plot(colSums(MG94.Force.nonclock.filtered[10:21, ]), MG94.nonclock.filtered[8, ],
     type = "n", xlab = "Total branch length", ylab = "Tau" )
col.color <- rep("black", dim(MG94.nonclock.filtered)[2])
col.color[(MG94.nonclock.filtered - MG94.Force.nonclock.filtered)[2,] < improvement.lmt] <- "red"
points(x = colSums(MG94.Force.nonclock.filtered[10:21, ]), y = MG94.nonclock.filtered[8, ],
       type = "p", col = col.color, bg = col.color)
legend("topright", 
       c(paste("lnL improvement > ", toString(improvement.lmt), sep = ""),
         paste("lnL improvement < ", toString(improvement.lmt), sep = "")),
       lty = c(1, 1), 
       lwd = c(2.5, 2.5),
       col = c("black", "red"))
title("MG94 nonclock")

MG94 clock case

improvement.lmt <- 20
plot(colSums(MG94.Force.clock.filtered[10:21, ]), MG94.clock.filtered[8, ],
     type = "n", xlab = "Total branch length", ylab = "Tau" )
col.color <- rep("black", dim(MG94.clock.filtered)[2])
col.color[(MG94.clock.filtered - MG94.Force.clock.filtered)[2,] < improvement.lmt] <- "red"
points(x = colSums(MG94.Force.clock.filtered[10:21, ]), y = MG94.clock.filtered[8, ],
       type = "p", col = col.color, bg = col.color)
legend("topright", 
       c(paste("lnL improvement > ", toString(improvement.lmt), sep = ""),
         paste("lnL improvement < ", toString(improvement.lmt), sep = "")),
       lty = c(1, 1), 
       lwd = c(2.5, 2.5),
       col = c("black", "red"))
title("MG94 clock")

Now see if the pairs red in HKY are also red in MG94

nonclock case

HKY.pair.names[(HKY.nonclock.filtered - HKY.Force.nonclock.filtered)[2,] < 50]

##  [1] "YLR406C_YDL075W" "YIL057C_YER067W" "YPL087W_YBR183W"
##  [4] "YDR518W_YCL043C" "YBR024W_YBR037C" "YPL037C_YDR252W"
##  [7] "YGL133W_YPL216W" "YLR284C_YOR180C" "YAL056W_YOR371C"
## [10] "YDR438W_YML018C"

MG94.pair.names[(MG94.nonclock.filtered - MG94.Force.nonclock.filtered)[2,] < improvement.lmt]

##  [1] "YLR406C_YDL075W" "YER131W_YGL189C" "YDR502C_YLR180W"
##  [4] "YIL057C_YER067W" "YPL087W_YBR183W" "YGR043C_YLR354C"
##  [7] "YPR157W_YGR141W" "YDR099W_YER177W" "YBR024W_YBR037C"
## [10] "YPL037C_YDR252W" "YGL133W_YPL216W" "YBR117C_YPR074C"
## [13] "YLR284C_YOR180C" "YMR243C_YOR316C" "YAL056W_YOR371C"
## [16] "YDR438W_YML018C"

Only 3 pairs show up in both: YLR406C_YDL075W, YIL057C_YER067W, YDR438W_YML018C.

clock case

HKY.pair.names[(HKY.clock.filtered - HKY.Force.clock.filtered)[2,] < 50]

##  [1] "YLR406C_YDL075W" "YER131W_YGL189C" "YIL057C_YER067W"
##  [4] "YPL087W_YBR183W" "YDR518W_YCL043C" "YGR043C_YLR354C"
##  [7] "YBR024W_YBR037C" "YPL037C_YDR252W" "YGL133W_YPL216W"
## [10] "YLR284C_YOR180C" "YAL056W_YOR371C" "YDR438W_YML018C"

MG94.pair.names[(MG94.clock.filtered - MG94.Force.clock.filtered)[2,] < improvement.lmt]

##  [1] "YLR406C_YDL075W" "YER131W_YGL189C" "YIL057C_YER067W"
##  [4] "YPL087W_YBR183W" "YGR043C_YLR354C" "YPR157W_YGR141W"
##  [7] "YDR099W_YER177W" "YBR024W_YBR037C" "YPL037C_YDR252W"
## [10] "YGL133W_YPL216W" "YBR117C_YPR074C" "YLR284C_YOR180C"
## [13] "YMR243C_YOR316C" "YAL056W_YOR371C" "YDR438W_YML018C"

=========================================

Now plot Total blen v.s. Tau into different groups (differ by lnL improvement per site)

HKY nonclock case

plot(colSums(HKY.Force.nonclock.filtered[9:20, ]), HKY.nonclock.filtered[8, ],
     type = "n", xlab = "Total branch length", ylab = "Tau" )
col.color <- rep("black", dim(HKY.nonclock.filtered)[2])
col.color[((HKY.nonclock.filtered - HKY.Force.nonclock.filtered)[2,] /  HKY.nonclock.filtered[1, ]) < 0.06] <- "red"
points(x = colSums(HKY.Force.nonclock.filtered[9:20, ]), y = HKY.nonclock.filtered[8, ],
       type = "p", col = col.color, bg = col.color)
legend("topright", c("lnL improvement per site > 0.06",
                     "lnL improvement per site < 0.06"),
       lty = c(1, 1), 
       lwd = c(2.5, 2.5),
       col = c("black", "red"))
title("HKY nonclock")

# lnL improvement per site
((HKY.nonclock.filtered - HKY.Force.nonclock.filtered)[2,] /  HKY.nonclock.filtered[1, ])

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##        0.134895        0.153594        0.185235        0.323550 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##        0.274350        0.173246        0.090389        0.050645 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##        0.152353        0.028734        0.060532        0.060868 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##        0.210301        0.053765        0.031368        0.234185 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##        0.115064        0.163738        0.009705        0.229431 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##        0.252429        0.389454        0.196469        0.063228 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##        0.047468        0.127221        0.020901        0.281591 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##        0.394727        0.037829        0.076393        0.338947 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##        0.344680        0.390797        0.011673        0.017508

HKY clock case

plot(colSums(HKY.Force.clock.filtered[9:20, ]), HKY.clock.filtered[8, ],
     type = "n", xlab = "Total branch length", ylab = "Tau" )
col.color <- rep("black", dim(HKY.clock.filtered)[2])
col.color[((HKY.clock.filtered - HKY.Force.clock.filtered)[2,] /  HKY.clock.filtered[1, ]) < 0.06] <- "red"
points(x = colSums(HKY.Force.clock.filtered[9:20, ]), y = HKY.clock.filtered[8, ],
       type = "p", col = col.color, bg = col.color)
legend("topright", c("lnL improvement per site > 0.06",
                     "lnL improvement per site < 0.06"),
       lty = c(1, 1), 
       lwd = c(2.5, 2.5),
       col = c("black", "red"))
title("HKY clock")

# lnL improvement per site
((HKY.clock.filtered - HKY.Force.clock.filtered)[2,] /  HKY.clock.filtered[1, ])

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##        0.135924        0.133484        0.148359        0.322104 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##        0.269324        0.152208        0.088402        0.043125 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##        0.131456        0.022857        0.040853        0.052711 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##        0.204649        0.038632        0.035687        0.234133 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##        0.115533        0.167633        0.004019        0.229508 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##        0.249856        0.410427        0.191646        0.048577 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##        0.037242        0.106903        0.022112        0.278139 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##        0.372669        0.030555        0.059616        0.370227 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##        0.341618        0.385737        0.013551        0.007138

MG94 nonclock case

improvement.lmt <- 0.16
plot(colSums(MG94.Force.nonclock.filtered[10:21, ]), MG94.nonclock.filtered[8, ],
     type = "n", xlab = "Total branch length", ylab = "Tau" )
col.color <- rep("black", dim(MG94.nonclock.filtered)[2])
col.color[((MG94.nonclock.filtered - MG94.Force.nonclock.filtered)[2,] / MG94.nonclock.filtered[1]) < improvement.lmt] <- "red"
points(x = colSums(MG94.Force.nonclock.filtered[10:21, ]), y = MG94.nonclock.filtered[8, ],
       type = "p", col = col.color, bg = col.color)
legend("topright", 
       c(paste("lnL improvement per site > ", toString(improvement.lmt), sep = ""),
         paste("lnL improvement per site < ", toString(improvement.lmt), sep = "")),
       lty = c(1, 1), 
       lwd = c(2.5, 2.5),
       col = c("black", "red"))
title("MG94 nonclock")

# lnL improvement per site
((MG94.nonclock.filtered - MG94.Force.nonclock.filtered)[2,] / MG94.nonclock.filtered[1])

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##        0.151632        0.060301        0.050936        0.504254 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##        0.676453        0.198674        0.146466        0.105363 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##        0.523279        0.301559        0.049232        0.027064 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##        0.127278        0.072056        0.033617        0.332390 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##        0.486996        0.326961        0.048522        0.192824 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##        0.508448        0.978479        0.286073        0.469334 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##        0.163692        0.268073        0.007964        0.296685 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##        2.519933        0.276193        0.067156        0.746849 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##        1.242955        1.166951        0.130769        0.028611

MG94 clock case

improvement.lmt <- 0.4
plot(colSums(MG94.Force.clock.filtered[10:21, ]), MG94.clock.filtered[8, ],
     type = "n", xlab = "Total branch length", ylab = "Tau" )
col.color <- rep("black", dim(MG94.clock.filtered)[2])
col.color[((MG94.clock.filtered - MG94.Force.clock.filtered)[2,] / MG94.clock.filtered[1]) < improvement.lmt] <- "red"
points(x = colSums(MG94.Force.clock.filtered[10:21, ]), y = MG94.clock.filtered[8, ],
       type = "p", col = col.color, bg = col.color)
legend("topright", 
       c(paste("lnL improvement per site > ", toString(improvement.lmt), sep = ""),
         paste("lnL improvement per site < ", toString(improvement.lmt), sep = "")),
       lty = c(1, 1), 
       lwd = c(2.5, 2.5),
       col = c("black", "red"))
title("MG94 clock")

# lnL improvement per site
((MG94.clock.filtered - MG94.Force.clock.filtered)[2,] / MG94.clock.filtered[1])

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##        0.155574        0.138675        0.352552        0.609461 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##        0.709713        0.181451        0.149269        0.101635 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##        0.475907        0.274049        0.028466        0.030841 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##        0.140684        0.059552        0.054233        0.348601 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##        0.514224        0.345646        0.016160        0.209024 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##        0.597249        1.053827        0.272536        0.429067 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##        0.121029        0.257700        0.007321        0.442413 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##        4.286240        0.226330        0.048014        0.824266 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##        1.307444        1.258567        0.162305        0.005364

OK, the MG94 nonclock long branch lengths are suspecious. It seems that the first four branches carry most of the total length. What if plot without those branches?

MG94 nonclock case without first four branches: (N0, N1), (N0,kluyveri), (N1,N2), (N1,castellii) are excluded

improvement.lmt <- 0.16
plot(colSums(MG94.Force.nonclock.filtered[14:21, ]), MG94.nonclock.filtered[8, ],
     type = "n", xlab = "Partial branch length", ylab = "Tau" )
col.color <- rep("black", dim(MG94.nonclock.filtered)[2])
col.color[((MG94.nonclock.filtered - MG94.Force.nonclock.filtered)[2,] / MG94.nonclock.filtered[1]) < improvement.lmt] <- "red"
points(x = colSums(MG94.Force.nonclock.filtered[14:21, ]), y = MG94.nonclock.filtered[8, ],
       type = "p", col = col.color, bg = col.color)
legend("topright", 
       c(paste("lnL improvement per site > ", toString(improvement.lmt), sep = ""),
         paste("lnL improvement per site < ", toString(improvement.lmt), sep = "")),
       lty = c(1, 1), 
       lwd = c(2.5, 2.5),
       col = c("black", "red"))
title("MG94 nonclock 4 branches removed")

# lnL improvement per site
((MG94.nonclock.filtered - MG94.Force.nonclock.filtered)[2,] / MG94.nonclock.filtered[1])

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##        0.151632        0.060301        0.050936        0.504254 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##        0.676453        0.198674        0.146466        0.105363 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##        0.523279        0.301559        0.049232        0.027064 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##        0.127278        0.072056        0.033617        0.332390 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##        0.486996        0.326961        0.048522        0.192824 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##        0.508448        0.978479        0.286073        0.469334 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##        0.163692        0.268073        0.007964        0.296685 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##        2.519933        0.276193        0.067156        0.746849 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##        1.242955        1.166951        0.130769        0.028611

=========================================

Now, plot posterior expected Number of geneconv events from paralog 1 to paralog2 / paralog2 to paralog 1 of each pair

Each plot is for each pair of paralogs under one model with each point representing one branch

# First read in descriptions of the cerevisiae genes.
description <- as.list(readLines("/Users/xji3/Genconv/NewClusterPackRun/cerevisiae_genes_description.txt"))
names(description) <- readLines("/Users/xji3/Genconv/NewClusterPackRun/cerevisiae_genes.txt")

for (i in 1:36){

#HKY nonclock

  rows <- 35:44
  paralog <- strsplit(colnames(HKY.nonclock.filtered)[i], split = "_")[[1]]
  plot(HKY.nonclock.filtered[rows, i],  HKY.nonclock.filtered[rows + 12, i],
       xlab = paste(paralog[1], paralog[2], sep = " -> "),
       ylab = paste(paralog[2], paralog[1], sep = " -> "),
       main = paste(colnames(HKY.nonclock.filtered)[i], " HKY Non clock"))
  abline(a = 0, b = 1)


#HKY clock

  paralog <- strsplit(colnames(HKY.clock.filtered)[i], split = "_")[[1]]
  plot(HKY.clock.filtered[rows, i],  HKY.clock.filtered[rows + 12, i],
       xlab = paste(paralog[1], paralog[2], sep = " -> "),
       ylab = paste(paralog[2], paralog[1], sep = " -> "),
       main = paste(colnames(HKY.clock.filtered)[i], " HKY Clock"))
  abline(a = 0, b = 1)

#MG94 Non-clock

  rows <- 36:45
  paralog <- strsplit(colnames(MG94.nonclock.filtered)[i], split = "_")[[1]]
  plot(MG94.nonclock.filtered[rows, i],  MG94.nonclock.filtered[rows + 12, i],
       xlab = paste(paralog[1], paralog[2], sep = " -> "),
       ylab = paste(paralog[2], paralog[1], sep = " -> "),
       main = paste(colnames(MG94.nonclock.filtered)[i], " MG94 Non clock"))
  abline(a = 0, b = 1)

#MG94 clock

  paralog <- strsplit(colnames(MG94.clock.filtered)[i], split = "_")[[1]]
  plot(MG94.clock.filtered[rows, i],  MG94.clock.filtered[rows + 12, i],
       xlab = paste(paralog[1], paralog[2], sep = " -> "),
       ylab = paste(paralog[2], paralog[1], sep = " -> "),
       main = paste(colnames(MG94.clock.filtered)[i], " MG94 Clock"))
  abline(a = 0, b = 1)

# Now show the description from cerevisiae genes of this pair
print(description[paralog[1]][[1]])
print(description[paralog[2]][[1]])
}

## [1] "YLR406C RPL31B SGDID:S000004398, Chr XII from 931349-931065,931755-931699, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl31Ap and has similarity to rat L31 ribosomal protein; associates with the karyopherin Sxm1p; loss of both Rpl31p and Rpl39p confers lethality\""
## [1] "YDL075W RPL31A SGDID:S000002233, Chr IV from 322226-322282,322704-322988, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl31Bp and has similarity to rat L31 ribosomal protein; associates with the karyopherin Sxm1p; loss of both Rpl31p and Rpl39p confers lethality\""

## [1] "YER131W RPS26B SGDID:S000000933, Chr V from 423952-424311, Genome Release 64-1-1, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; nearly identical to Rps26Ap and has similarity to rat S26 ribosomal protein\""
## [1] "YGL189C RPS26A SGDID:S000003157, Chr VII from 148588-148229, Genome Release 64-1-1, reverse complement, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; nearly identical to Rps26Bp and has similarity to rat S26 ribosomal protein\""

## [1] "YDR502C SAM2 SGDID:S000002910, Chr IV from 1454464-1453310, Genome Release 64-1-1, reverse complement, Verified ORF, \"S-adenosylmethionine synthetase, catalyzes transfer of the adenosyl group of ATP to the sulfur atom of methionine; one of two differentially regulated isozymes (Sam1p and Sam2p)\""
## [1] "YLR180W SAM1 SGDID:S000004170, Chr XII from 515262-516410, Genome Release 64-1-1, Verified ORF, \"S-adenosylmethionine synthetase, catalyzes transfer of the adenosyl group of ATP to the sulfur atom of methionine; one of two differentially regulated isozymes (Sam1p and Sam2p)\""

## [1] "YML026C RPS18B SGDID:S000004488, Chr XIII from 223380-222987,223828-223782, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; nearly identical to Rps18Ap and has similarity to E. coli S13 and rat S18 ribosomal proteins\""
## [1] "YDR450W RPS18A SGDID:S000002858, Chr IV from 1359923-1359969,1360405-1360798, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; nearly identical to Rps18Bp and has similarity to E. coli S13 and rat S18 ribosomal proteins\""

## [1] "YNL301C RPL18B SGDID:S000005245, Chr XIV from 64018-63570,64562-64451, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, identical to Rpl18Ap and has similarity to rat L18 ribosomal protein\""
## [1] "YOL120C RPL18A SGDID:S000005480, Chr XV from 93843-93395,94402-94291, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, identical to Rpl18Bp and has similarity to rat L18 ribosomal protein; intron of RPL18A pre-mRNA forms stem-loop structures that are a target for Rnt1p cleavage leading to degradation\""

## [1] "YHR106W TRR2 SGDID:S000001148, Chr VIII from 325598-326626, Genome Release 64-1-1, Verified ORF, \"Mitochondrial thioredoxin reductase involved in protection against oxidative stress, required with Glr1p to maintain the redox state of Trx3p; contains active-site motif (CAVC) present in prokaryotic orthologs; binds NADPH and FAD\""
## [1] "YDR353W TRR1 SGDID:S000002761, Chr IV from 1183299-1184258, Genome Release 64-1-1, Verified ORF, \"Cytoplasmic thioredoxin reductase, key regulatory enzyme that determines the redox state of the thioredoxin system, which acts as a disulfide reductase system and protects cells against both oxidative and reductive stress\""

## [1] "YIL057C RGI2 SGDID:S000001319, Chr IX from 248396-247902, Genome Release 64-1-1, reverse complement, Verified ORF, \"Protein of unknown function involved in energy metabolism under respiratory conditions; expression induced under carbon limitation and repressed under high glucose\""
## [1] "YER067W RGI1 SGDID:S000000869, Chr V from 292066-292551, Genome Release 64-1-1, Verified ORF, \"Protein of unknown function involved in energy metabolism under respiratory conditions; protein abundance is increased upon intracellular iron depletion\""

## [1] "YPL087W YDC1 SGDID:S000006008, Chr XVI from 383455-384408, Genome Release 64-1-1, Verified ORF, \"Alkaline dihydroceramidase, involved in sphingolipid metabolism; preferentially hydrolyzes dihydroceramide to a free fatty acid and dihydrosphingosine; has a minor reverse activity\""
## [1] "YBR183W YPC1 SGDID:S000000387, Chr II from 596115-597065, Genome Release 64-1-1, Verified ORF, \"Alkaline ceramidase that also has reverse (CoA-independent) ceramide synthase activity, catalyzes both breakdown and synthesis of phytoceramide; overexpression confers fumonisin B1 resistance\""

## [1] "YNL069C RPL16B SGDID:S000005013, Chr XIV from 494524-493956,495001-494974, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"N-terminally acetylated protein component of the large (60S) ribosomal subunit, binds to 5.8 S rRNA; has similarity to Rpl16Ap, E. coli L13 and rat L13a ribosomal proteins; transcriptionally regulated by Rap1p\""
## [1] "YIL133C RPL16A SGDID:S000001395, Chr IX from 99095-98527,99416-99386, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"N-terminally acetylated protein component of the large (60S) ribosomal subunit, binds to 5.8 S rRNA; has similarity to Rpl16Bp, E. coli L13 and rat L13a ribosomal proteins; transcriptionally regulated by Rap1p\""

## [1] "YDR518W EUG1 SGDID:S000002926, Chr IV from 1478608-1480161, Genome Release 64-1-1, Verified ORF, \"Protein disulfide isomerase of the endoplasmic reticulum lumen, function overlaps with that of Pdi1p; may interact with nascent polypeptides in the ER\""
## [1] "YCL043C PDI1 SGDID:S000000548, Chr III from 50221-48653, Genome Release 64-1-1, reverse complement, Verified ORF, \"Protein disulfide isomerase; multifunctional protein resident in the endoplasmic reticulum lumen, essential for the formation of disulfide bonds in secretory and cell-surface proteins, unscrambles non-native disulfide bonds; forms a complex with Mnl1p that has exomannosidase activity, processing unfolded protein-bound Man8GlcNAc2 oligosaccharides to Man7GlcNAc2 which promotes degradation in the unfolded protein response\""

## [1] "YGR043C NQM1 SGDID:S000003275, Chr VII from 581436-580435, Genome Release 64-1-1, reverse complement, Verified ORF, \"Transaldolase of unknown function; transcription is repressed by Mot1p and induced by alpha-factor and during diauxic shift\""
## [1] "YLR354C TAL1 SGDID:S000004346, Chr XII from 837357-836350, Genome Release 64-1-1, reverse complement, Verified ORF, \"Transaldolase, enzyme in the non-oxidative pentose phosphate pathway; converts sedoheptulose 7-phosphate and glyceraldehyde 3-phosphate to erythrose 4-phosphate and fructose 6-phosphate\""

## [1] "YPR157W TDA6 SGDID:S000006361, Chr XVI from 841266-842669, Genome Release 64-1-1, Uncharacterized ORF, \"Putative protein of unknown function; induced by treatment with 8-methoxypsoralen and UVA irradiation; null mutant is sensitive to expression of the top1-T722A allele\""
## [1] "YGR141W VPS62 SGDID:S000003373, Chr VII from 770569-771972, Genome Release 64-1-1, Verified ORF, \"Vacuolar protein sorting (VPS) protein required for cytoplasm to vacuole targeting of proteins\""

## [1] "YDR099W BMH2 SGDID:S000002506, Chr IV from 653607-654428, Genome Release 64-1-1, Verified ORF, \"14-3-3 protein, minor isoform; controls proteome at post-transcriptional level, binds proteins and DNA, involved in regulation of many processes including exocytosis, vesicle transport, Ras/MAPK signaling, and rapamycin-sensitive signaling\""
## [1] "YER177W BMH1 SGDID:S000000979, Chr V from 545611-546414, Genome Release 64-1-1, Verified ORF, \"14-3-3 protein, major isoform; controls proteome at post-transcriptional level, binds proteins and DNA, involved in regulation of many processes including exocytosis, vesicle transport, Ras/MAPK signaling, and rapamycin-sensitive signaling\""

## [1] "YBR024W SCO2 SGDID:S000000228, Chr II from 289445-290350, Genome Release 64-1-1, Verified ORF, \"Protein anchored to the mitochondrial inner membrane, similar to Sco1p and may have a redundant function with Sco1p in delivery of copper to cytochrome c oxidase; interacts with Cox2p\""
## [1] "YBR037C SCO1 SGDID:S000000241, Chr II from 311452-310565, Genome Release 64-1-1, reverse complement, Verified ORF, \"Copper-binding protein of the mitochondrial inner membrane, required for cytochrome c oxidase activity and respiration; may function to deliver copper to cytochrome c oxidase; has similarity to thioredoxins\""

## [1] "YPL037C EGD1 SGDID:S000005958, Chr XVI from 481901-481428, Genome Release 64-1-1, reverse complement, Verified ORF, \"Subunit beta1 of the nascent polypeptide-associated complex (NAC) involved in protein targeting, associated with cytoplasmic ribosomes; enhances DNA binding of the Gal4p activator; homolog of human BTF3b\""
## [1] "YDR252W BTT1 SGDID:S000002660, Chr IV from 963412-963861, Genome Release 64-1-1, Verified ORF, \"Beta3 subunit of the heterotrimeric nascent polypeptide-associated complex which binds ribosomes via its beta-subunits in close proximity to nascent polypeptides; interacts with Caf130p of the CCR4-NOT complex; similar to human BTF3\""

## [1] "YMR143W RPS16A SGDID:S000004751, Chr XIII from 551928-551951,552496-552903, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; identical to Rps16Bp and has similarity to E. coli S9 and rat S16 ribosomal proteins\""
## [1] "YDL083C RPS16B SGDID:S000002241, Chr IV from 307333-306926,307789-307766, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; identical to Rps16Ap and has similarity to E. coli S9 and rat S16 ribosomal proteins\""

## [1] "YPR159W KRE6 SGDID:S000006363, Chr XVI from 857583-859745, Genome Release 64-1-1, Verified ORF, \"Type II integral membrane protein required for beta-1,6 glucan biosynthesis; putative beta-glucan synthase; localizes to the ER, plasma membrane, sites of polarized growth and secretory vesicles; functionally redundant with Skn1p\""
## [1] "YGR143W SKN1 SGDID:S000003375, Chr VII from 775193-777508, Genome Release 64-1-1, Verified ORF, \"Protein involved in sphingolipid biosynthesis; type II membrane protein with similarity to Kre6p\""

## [1] "YJL177W RPL17B SGDID:S000003713, Chr X from 90786-91094,91412-91657, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl17Ap and has similarity to E. coli L22 and rat L17 ribosomal proteins\""
## [1] "YKL180W RPL17A SGDID:S000001663, Chr XI from 109269-109577,109884-110129, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl17Bp and has similarity to E. coli L22 and rat L17 ribosomal proteins; copurifies with the Dam1 complex (aka DASH complex)\""

## [1] "YGL133W ITC1 SGDID:S000003101, Chr VII from 257707-261501, Genome Release 64-1-1, Verified ORF, \"Subunit of the ATP-dependent Isw2p-Itc1p chromatin remodeling complex, required for repression of a-specific genes, repression of early meiotic genes during mitotic growth, and repression of INO1; similar to mammalian Acf1p, the regulatory subunit of the mammalian ATP-utilizing chromatin assembly and modifying factor (ACF) complex\""
## [1] "YPL216W YPL216W SGDID:S000006137, Chr XVI from 143821-147129, Genome Release 64-1-1, Uncharacterized ORF, \"Putative protein of unknown function; YPL216W is not an essential gene\""

## [1] "YJR048W CYC1 SGDID:S000003809, Chr X from 526335-526664, Genome Release 64-1-1, Verified ORF, \"Cytochrome c, isoform 1; electron carrier of the mitochondrial intermembrane space that transfers electrons from ubiquinone-cytochrome c oxidoreductase to cytochrome c oxidase during cellular respiration\""
## [1] "YEL039C CYC7 SGDID:S000000765, Chr V from 79977-79636, Genome Release 64-1-1, reverse complement, Verified ORF, \"Cytochrome c isoform 2, expressed under hypoxic conditions; electron carrier of the mitochondrial intermembrane space that transfers electrons from ubiquinone-cytochrome c oxidoreductase to cytochrome c oxidase during cellular respiration\""

## [1] "YBR191W RPL21A SGDID:S000000395, Chr II from 606270-606280,606669-607140, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl21Bp and has similarity to rat L21 ribosomal protein\""
## [1] "YPL079W RPL21B SGDID:S000006000, Chr XVI from 406636-406646,407068-407539, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl21Ap and has similarity to rat L21 ribosomal protein\""

## [1] "YER074W RPS24A SGDID:S000000876, Chr V from 306323-306325,306792-307196, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; identical to Rps24Bp and has similarity to rat S24 ribosomal protein\""
## [1] "YIL069C RPS24B SGDID:S000001331, Chr IX from 231957-231553,232369-232367, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; identical to Rps24Ap and has similarity to rat S24 ribosomal protein\""

## [1] "YDR418W RPL12B SGDID:S000002826, Chr IV from 1301616-1302113, Genome Release 64-1-1, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl12Ap; rpl12a rpl12b double mutant exhibits slow growth and slow translation; has similarity to E. coli L11 and rat L12 ribosomal proteins\""
## [1] "YEL054C RPL12A SGDID:S000000780, Chr V from 53218-52721, Genome Release 64-1-1, reverse complement, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl12Bp; rpl12a rpl12b double mutant exhibits slow growth and slow translation; has similarity to E. coli L11 and rat L12 ribosomal proteins\""

## [1] "YNL049C SFB2 SGDID:S000004994, Chr XIV from 537911-535281, Genome Release 64-1-1, reverse complement, Verified ORF, \"Component of the Sec23p-Sfb2p heterodimer of the COPII vesicle coat, required for cargo selection during vesicle formation in ER to Golgi transport; homologous to Sec24p and Sfb3p\""
## [1] "YIL109C SEC24 SGDID:S000001371, Chr IX from 160165-157385, Genome Release 64-1-1, reverse complement, Verified ORF, \"Component of the Sec23p-Sec24p heterodimer of the COPII vesicle coat, required for cargo selection during vesicle formation in ER to Golgi transport; homologous to Sfb2p and Sfb3p\""

## [1] "YBR117C TKL2 SGDID:S000000321, Chr II from 476437-474392, Genome Release 64-1-1, reverse complement, Verified ORF, \"Transketolase, similar to Tkl1p; catalyzes conversion of xylulose-5-phosphate and ribose-5-phosphate to sedoheptulose-7-phosphate and glyceraldehyde-3-phosphate in the pentose phosphate pathway; needed for synthesis of aromatic amino acids\""
## [1] "YPR074C TKL1 SGDID:S000006278, Chr XVI from 694838-692796, Genome Release 64-1-1, reverse complement, Verified ORF, \"Transketolase, similar to Tkl2p; catalyzes conversion of xylulose-5-phosphate and ribose-5-phosphate to sedoheptulose-7-phosphate and glyceraldehyde-3-phosphate in the pentose phosphate pathway; needed for synthesis of aromatic amino acids\""

## [1] "YPL232W SSO1 SGDID:S000006153, Chr XVI from 107275-108147, Genome Release 64-1-1, Verified ORF, \"Plasma membrane t-SNARE involved in fusion of secretory vesicles at the plasma membrane and in vesicle fusion during sporulation; forms a complex with Sec9p that binds v-SNARE Snc2p; syntaxin homolog; functionally redundant with Sso2p\""
## [1] "YMR183C SSO2 SGDID:S000004795, Chr XIII from 627808-626921, Genome Release 64-1-1, reverse complement, Verified ORF, \"Plasma membrane t-SNARE involved in fusion of secretory vesicles at the plasma membrane; syntaxin homolog that is functionally redundant with Sso1p\""

## [1] "YLR284C ECI1 SGDID:S000004274, Chr XII from 707040-706198, Genome Release 64-1-1, reverse complement, Verified ORF, \"Peroxisomal delta3,delta2-enoyl-CoA isomerase, hexameric protein that converts 3-hexenoyl-CoA to trans-2-hexenoyl-CoA, essential for the beta-oxidation of unsaturated fatty acids, oleate-induced\""
## [1] "YOR180C DCI1 SGDID:S000005706, Chr XV from 675167-674352, Genome Release 64-1-1, reverse complement, Verified ORF, \"Peroxisomal protein; identification as a delta(3,5)-delta(2,4)-dienoyl-CoA isomerase involved in fatty acid metabolism is disputed\""

## [1] "YBL087C RPL23A SGDID:S000000183, Chr II from 60193-59822,60739-60698, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, identical to Rpl23Bp and has similarity to E. coli L14 and rat L23 ribosomal proteins\""
## [1] "YER117W RPL23B SGDID:S000000919, Chr V from 396769-396810,397282-397653, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, identical to Rpl23Ap and has similarity to E. coli L14 and rat L23 ribosomal proteins\""

## [1] "YGL062W PYC1 SGDID:S000003030, Chr VII from 385196-388732, Genome Release 64-1-1, Verified ORF, \"Pyruvate carboxylase isoform, cytoplasmic enzyme that converts pyruvate to oxaloacetate; highly similar to isoform Pyc2p but differentially regulated; mutations in the human homolog are associated with lactic acidosis\""
## [1] "YBR218C PYC2 SGDID:S000000422, Chr II from 662249-658707, Genome Release 64-1-1, reverse complement, Verified ORF, \"Pyruvate carboxylase isoform, cytoplasmic enzyme that converts pyruvate to oxaloacetate; highly similar to isoform Pyc1p but differentially regulated; mutations in the human homolog are associated with lactic acidosis\""

## [1] "YIR033W MGA2 SGDID:S000001472, Chr IX from 416124-419465, Genome Release 64-1-1, Verified ORF, \"ER membrane protein involved in regulation of OLE1 transcription, acts with homolog Spt23p; inactive ER form dimerizes and one subunit is then activated by ubiquitin/proteasome-dependent processing followed by nuclear targeting\""
## [1] "YKL020C SPT23 SGDID:S000001503, Chr XI from 402080-398832, Genome Release 64-1-1, reverse complement, Verified ORF, \"ER membrane protein involved in regulation of OLE1 transcription, acts with homolog Mga2p; inactive ER form dimerizes and one subunit is then activated by ubiquitin/proteasome-dependent processing followed by nuclear targeting\""

## [1] "YMR243C ZRC1 SGDID:S000004856, Chr XIII from 756166-754838, Genome Release 64-1-1, reverse complement, Verified ORF, \"Vacuolar membrane zinc transporter, transports zinc from the cytosol into the vacuole for storage; also has a role in resistance to zinc shock resulting from a sudden influx of zinc into the cytoplasm\""
## [1] "YOR316C COT1 SGDID:S000005843, Chr XV from 907555-906236, Genome Release 64-1-1, reverse complement, Verified ORF, \"Vacuolar transporter that mediates zinc transport into the vacuole; overexpression confers resistance to cobalt and rhodium\""

## [1] "YLR333C RPS25B SGDID:S000004325, Chr XII from 795899-795573, Genome Release 64-1-1, reverse complement, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; nearly identical to Rps25Ap and has similarity to rat S25 ribosomal protein\""
## [1] "YGR027C RPS25A SGDID:S000003259, Chr VII from 534458-534132, Genome Release 64-1-1, reverse complement, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; nearly identical to Rps25Bp and has similarity to rat S25 ribosomal protein\""

## [1] "YMR142C RPL13B SGDID:S000004750, Chr XIII from 550801-550206,551207-551204, Genome Release 64-1-1, reverse complement, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl13Ap; not essential for viability; has similarity to rat L13 ribosomal protein\""
## [1] "YDL082W RPL13A SGDID:S000002240, Chr IV from 308424-308427,308793-309388, Genome Release 64-1-1, intron sequence removed, Verified ORF, \"Protein component of the large (60S) ribosomal subunit, nearly identical to Rpl13Bp; not essential for viability; has similarity to rat L13 ribosomal protein\""

## [1] "YER102W RPS8B SGDID:S000000904, Chr V from 363100-363702, Genome Release 64-1-1, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; identical to Rps8Ap and has similarity to rat S8 ribosomal protein\""
## [1] "YBL072C RPS8A SGDID:S000000168, Chr II from 89125-88523, Genome Release 64-1-1, reverse complement, Verified ORF, \"Protein component of the small (40S) ribosomal subunit; identical to Rps8Bp and has similarity to rat S8 ribosomal protein\""

## [1] "YAL056W GPB2 SGDID:S000000052, Chr I from 39259-41901, Genome Release 64-1-1, Verified ORF, \"Multistep regulator of cAMP-PKA signaling; inhibits PKA downstream of Gpa2p and Cyr1p, thereby increasing cAMP dependency; inhibits Ras activity through direct interactions with Ira1p/2p; regulated by G-alpha protein Gpa2p; homolog of Gpb1p\""
## [1] "YOR371C GPB1 SGDID:S000005898, Chr XV from 1034182-1031489, Genome Release 64-1-1, reverse complement, Verified ORF, \"Multistep regulator of cAMP-PKA signaling; inhibits PKA downstream of Gpa2p and Cyr1p, thereby increasing cAMP dependency; promotes ubiquitin-dependent proteolysis of Ira2p; regulated by G-alpha protein Gpa2p; homolog of Gpb2p\""

## [1] "YDR438W THI74 SGDID:S000002846, Chr IV from 1338274-1339386, Genome Release 64-1-1, Verified ORF, \"Mitochondrial transporter repressible by thiamine\""
## [1] "YML018C YML018C SGDID:S000004480, Chr XIII from 235952-234771, Genome Release 64-1-1, reverse complement, Uncharacterized ORF, \"Putative protein of unknown function; green fluorescent protein (GFP)-fusion protein localizes to the membrane of the vacuole; physical interaction with Atg27p suggests a possible role in autophagy; YML018C is not an essential gene\""

=========================================

Now, plot histogram of posterior expected Number of geneconv events from paralog 1 to paralog2 / paralog2 to paralog 1 of each pair

Each plot is for one model. Each point contains the sum of expected number of directional IGC over all branches.

# HKY nonclock
rows <- 35:44
plot(rowSums(HKY.nonclock.filtered[rows, 1:36]),
     rowSums(HKY.nonclock.filtered[rows + 12, 1:36]),
       xlab = paste("paralog 1", "paralog 2", sep = " -> "),
       ylab = paste("paralog 2", "paralog 1", sep = " -> "),
       main = paste("All pairs HKY nonclock"))
abline(a = 0, b = 1)

# HKY clock
plot(rowSums(HKY.clock.filtered[rows, 1:36]),
     rowSums(HKY.clock.filtered[rows + 12, 1:36]),
       xlab = paste("paralog 1", "paralog 2", sep = " -> "),
       ylab = paste("paralog 2", "paralog 1", sep = " -> "),
       main = paste("All pairs HKY clock"))
abline(a = 0, b = 1)

# MG94 nonclock
rows <- 36:45
plot(rowSums(MG94.nonclock.filtered[rows, 1:36]),
     rowSums(MG94.nonclock.filtered[rows + 12, 1:36]),
       xlab = paste("paralog 1", "paralog 2", sep = " -> "),
       ylab = paste("paralog 2", "paralog 1", sep = " -> "),
       main = paste("All pairs MG94 nonclock"))
abline(a = 0, b = 1)

# MG94 clock
plot(rowSums(MG94.clock.filtered[rows, 1:36]),
     rowSums(MG94.clock.filtered[rows + 12, 1:36]),
       xlab = paste("paralog 1", "paralog 2", sep = " -> "),
       ylab = paste("paralog 2", "paralog 1", sep = " -> "),
       main = paste("All pairs MG94 clock"))
abline(a = 0, b = 1)

=========================================

Now, calculate the LR test statistics of adding one more directional tau parameter

# HKY nonclock
Tlr.HKY.nonclock <- 2 * (Dir_HKY_nonclock_summary[2, ] - HKY_nonclock_summary[2, ])
# HKY clock
Tlr.HKY.clock <- 2 * (Dir_HKY_clock_summary[2, ] - HKY_clock_summary[2, ])

# MG94 nonclock
Tlr.MG94.nonclock <- 2 * (Dir_MG94_nonclock_summary[2, ] - MG94_nonclock_summary[2, ])
# MG94 clock
Tlr.MG94.clock <- 2 * (Dir_MG94_clock_summary[2, ] - MG94_clock_summary[2, ])

Tlr.HKY.nonclock

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##       1.211e+01       1.386e+01       5.128e+00       1.152e-03 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##       4.792e+01       4.929e+00       1.798e-01       1.175e+00 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##       4.783e-02       5.049e+00       1.823e+01       4.510e+00 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##       1.894e+00       1.580e-02       2.685e+01       1.293e-02 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##       5.787e-01       2.456e+00       1.163e+02       1.489e-01 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##       2.082e-01       8.360e+00       1.147e-02       3.083e+01 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##       1.967e+01       7.606e-01       1.283e+00       4.424e+00 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##       2.917e+00       2.565e+00       8.183e+00       2.541e+00 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##       6.383e+00       1.421e+00       3.878e+00       1.872e+01

Tlr.HKY.clock

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##       1.220e+01       1.096e+01       3.412e+00       1.985e-02 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##       4.211e+01       4.647e+00       1.445e-01       2.608e+00 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##       6.643e-04       6.918e+00       1.191e+01       3.582e+00 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##       1.392e+00       7.523e-02       3.068e+01       1.242e-02 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##       3.734e-01       1.658e+00       8.522e+01       1.016e+00 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##       2.330e-02       8.031e+00       8.099e-02       2.441e+01 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##       1.667e+01       5.491e-01       7.247e-01       5.024e+00 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##       2.855e+00       1.859e+00       9.011e+00       3.392e+00 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##       5.431e+00       2.057e+00       3.948e+00       1.350e+01

Tlr.MG94.nonclock

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##       1.079e+01       9.323e+00       6.548e-01       1.342e-01 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##       1.325e+01       2.634e+00       9.410e-01       2.234e+00 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##      -2.025e+01       1.963e+00       8.470e+00       9.619e-01 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##      -3.721e-04       6.883e-02       1.484e+01       1.127e-01 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##       1.973e+00       4.666e-01       2.354e+01       6.339e-02 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##       4.489e-01      -1.156e+01       4.281e-01       2.186e+01 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##       2.063e+01       9.938e-01       8.249e-02       6.700e+00 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##       3.453e+00       9.447e+00       4.034e+00      -1.074e+01 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##       6.585e+00       1.445e+00       2.693e+00       5.490e+00

Tlr.MG94.clock

## YLR406C_YDL075W YER131W_YGL189C YDR502C_YLR180W YML026C_YDR450W 
##       10.045372        6.564622        0.737165        0.129719 
## YNL301C_YOL120C YHR106W_YDR353W YIL057C_YER067W YPL087W_YBR183W 
##       33.275759        2.622489        0.536157        2.712693 
## YNL069C_YIL133C YDR518W_YCL043C YGR043C_YLR354C YPR157W_YGR141W 
##        0.004496        3.386383        7.402037        0.212040 
## YDR099W_YER177W YBR024W_YBR037C YPL037C_YDR252W YMR143W_YDL083C 
##        0.007142        0.066007       16.141948        0.103398 
## YPR159W_YGR143W YJL177W_YKL180W YGL133W_YPL216W YJR048W_YEL039C 
##        1.424472        1.506427       16.166363        0.322980 
## YBR191W_YPL079W YER074W_YIL069C YDR418W_YEL054C YNL049C_YIL109C 
##        0.018798        7.145638        0.554162       22.913581 
## YBR117C_YPR074C YPL232W_YMR183C YLR284C_YOR180C YBL087C_YER117W 
##       20.217751        1.243511        0.052682        4.903412 
## YGL062W_YBR218C YIR033W_YKL020C YMR243C_YOR316C YLR333C_YGR027C 
##        5.328474        7.618701        3.896956        2.318428 
## YMR142C_YDL082W YER102W_YBL072C YAL056W_YOR371C YDR438W_YML018C 
##        4.782092        3.133447        1.878884        2.227691